Autonomous Motility of Active Filaments due to Spontaneous Flow-Symmetry 
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We simulate the nonlocal Stokesian hydrodynamics of an elastic filament which is active due to a 
permanent distribution of stresslets along its contour. A bending instability of an initially straight 
filament spontaneously breaks fiow symmetry and leads to autonomous filament motion which, 
depending on conformational symmetry can be translational or rotational. At high ratios of activity 
to elasticity, the linear instability develops into nonlinear fluctuating states with large amplitude 
deformations. The dynamics of these states can be qualitatively understood as a superposition 
of translational and rotational motion associated with fllament conformational modes of opposite 
symmetry. Our results can be tested in molecular-motor filament mixtures, synthetic chains of 
autocatalytic particles or other linearly connected systems where chemical energy is converted to 
mechanical energy in a fiuid environment. 
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Components which convert chemical energy to me- 
chanical energy internally are ubiquitous in biology. 
Common examples where this conversion leads to au- 
tonomous propulsion are molecular motors (at the sub- 
cellular level) and bacteria (at the cellular level) [1]. Re- 
cently, biomimetic elements which convert chemical en- 
ergy into translational [2] or rotational [3] motion have 
been realized in the laboratory. While the detailed mech- 
anisms leading to autonomous propulsion in these biolog- 
ical and soft matter systems show a wonderful variety [4] , 
their collective behavior tends to be universal and can be 
understood by appealing to symmetries and conservation 
laws [5] . This realization has led to many studies of the 
collective properties of suspensions of hydrodynamically 
interacting autonomously motile particles [6]. 

There are ample instances in biology, however, where 
the conversion of chemical to mechanical energy is not 
confined to a particle-like element but is, instead, dis- 
tributed over a line-like element. Such a situation arises, 
for example, in a microtubule with a row of molecular mo- 
tors converting energy while walking on it. The mechani- 
cal energy thus obtained not only produces motion of the 
motors but also generates reaction forces on the micro- 
tubule, which can deform elastically in response. Hydro- 
dynamic interactions between the motors and between 
segments of the microtubule must be taken into account 
since both are surrounded by a fluid. This combination 
of elasticity, autonomous motility through energy con- 
version and hydrodynamics is found in biomimetic con- 
texts as well. A recent example is provided by mixtures 
of motors which crosslink and walk on polymer bundles. 
A remarkable cilia-like beating phenomenon is observed 
in these systems [7] . A polymer in which the monomeric 
units are autocatalytic nanorods provides a nonbiological 
example of energy conversion on linear elastic elements. 



Though such elements are yet to be realized in the lab- 
oratory, active elements coupled to passive components 
through covalent bonds have been synthesized [2] and 
may lead to new kinds of nanomachincs [3] . 

Motivated by these biological and biomimetic exam- 
ples, we study, in this Letter, a semiflexible elastic fila- 
ment immersed in a viscous fluid with energy convert- 
ing "active" elements distributed along its length. We 
present an equation of motion for the filament that incor- 
porates the effects of nonlinear elastic deformation, active 
processes and nonlocal Stokesian hydrodynamic interac- 
tions. We use the lattice Boltzmann (LB) method to 
numerically solve the active filament equation of motion. 
Our simulations show that a straight active filament is 
linearly unstable to transverse perturbations. Depending 
on the symmetry of the perturbation, the hydrodynamic 
fiows produce translational or rotational motion of the 
filament. Conformational symmetry provides a qualita- 
tive way of understanding the dynamics of the nonlinear 
fiuctuating states that develop at high ratios of activity 
to elasticity. We describe these results and our model in 
detail below. 

Model: Our model for the active filament consists of 
N beads, with coordinates r„, interacting through a po- 
tential given by 
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The two-body harmonic spring potential Ua,{h^) = 
\k{hm — 6o)^ penalizes departures of bm, the modulus 



of the bond vector 
librium value of 6o. 



bm = |rm - rm+i|, from its equi- 
The three-body bending potential 
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%(bm,bm+i) = k(1 — cosc^m) penalizes departures of 
the angle c/jm between consecutive bond vectors from its 
equilibrium value of zero. The rigidity parameter R is 
related to the bending rigidity as k = boii. The repul- 
sive Lennard- Jones potential {/lj vanishes if the distance 
between beads r„j„ = |r,„ — r„| exceeds ctlj- The 71-th 
bead experiences a force f„ = —dU/drn when the fil- 
ament stretches or bends from its equilibrium position. 
With the above choice of potential the connected beads 
approximate an inextensible, semiflexible, self-avoiding 
filament. 

Active nonequilibrium processes, such as those that 
convert chemical energy to mechanical energy, are inter- 
nal to the fluid and hence cannot add net momentum 
to it. Then, the integral of the force density on a sur- 
face enclosing the active element must vanish. This is 
ensured if the active force density is the divergence of 
a stress. Since the active processes cannot add angular 
momentum to the fluid, the stress must be symmetric 
[8] . The most dominant Stokesian singularity with these 
properties is the stresslet [9]. There is a remaining free- 
dom of the sign of the stresslet and its angle relative to 
the filament. Motivated by the tangential stresses ex- 
erted by motors walking on microtubules [7], we choose 
the stresslet to be extensile and oriented along the in- 
stantaneous tangent t„ to the filament, 

cr„ = cro(t„t„ - I/d) (2) 

where d is the spatial dimension and ao > sets the scale 
of the activity. The results of other choices of sign and 
orientation will be presented elsewhere. 

Elastic forces and active stresses produce velocities 
in the fluid. In the Stokesian regime, the velocity in 
a three-dimensional unbounded fluid at location r pro- 
duced by a force f at the origin is Wa(r) = Oap{r)fp 
where Oa^(r) = [dap + fafp)/ii:rir is the Oseen ten- 
sor, Cartesian directions are indicated by Greek indices, 
•q is the fluid shear viscosity and fa = fa/r. Simi- 
larly, the velocity at location r produced by a stresslet 
<T at the origin is Vaij) = Dap-y{Y)(Jp-y where Daji-yir) = 
i^faSp^ + ?>fa'fpf^)/?>'K-qr'^ [10]. In the presence of rigid 
or periodic boundaries the tensors O and D must be 
replaced by the appropriate Green's functions of Stokes 
flow that vanish at the boundaries or have the periodicity 
of the domain [10]. Similarly, two-dimensional Green's 
functions must be used when studying the motion of fil- 
aments in planar flows [9] . The velocity of the n-th bead 
is obtained by summing the force and activity contribu- 
tions from all beads, including itself, to the fluid velocity 
at its location. An isolated spherical bead with a force f 
acquires a velocity /z f where n is its mobility. By symme- 
try, an isolated spherical bead with a stresslet cr cannot 
acquire a velocity. This gives the following equation of 
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FIG. 1. (Color online) Spontaneous symmetry breaking due 
to a linear instability, (a) A straight filament produces a flow 
symmetric about both its axis and its midpoint. This flow 
can only produce an extensile motion of the filament, (b) 
The straight conformation, however, is finearly unstable to 
transverse perturbations. Here the perturbation is symmetric 
about an axis passing through the filament midpoint. This 
leads to a fiow which is no longer symmetric about the ini- 
tial axis. It enhances the instability and can produce center 
of mass motion of the filament. Bead positions are shown 
as yeUow filled circles, while the fluid velocity is shown as 
streamlines with the background color indicating its magni- 
tude. 

motion for the active filament: 

N 
m— 1 

where, for m = n, Oap = fJ-Sap and Dap-y = 0. Equa- 
tions (1), (2) and (3) represent our model for the nonlocal 
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FIG. 2. (Color) Filament motion for high activity, (a) For an even initial conformation with A — 250.7, the motion is 
predominantly translational. (b) For an identical initial conformation with A — 150.4, the motion is translational as well as 
rotational due to the spontaneous appearance of conformations with odd symmetry. In both (a) and (b) the color of the trace 
corresponds to the center of mass speed |R| normalized by its maximum value, (c) and (d) Time traces of the x (top) and 
y components (bottom) of R (olive dashed line) and K — — {ao / Anrjbo) (xrn) (black solid line), where (xrn) is the average 
curvature, corresponding to simulations in (a) and (b) respectively. The velocity and the curvature remain strongly correlated 
for various filament conformational modes and activity numbers. Times are in 10^ LB steps. 



Stokesian hydrodynamics of an active elastic filament. In 
the absence of bending rigidity and activity, our model 
reduces to Zimm dynamics of a polymer in a good solvent 
[11]. 

The ratio of the strcsslct and Stokeslet terms in the 
equation of motion is a dimcnsionless measure of activity. 
Estimating the curvature elastic force as k/L^, where 
L = {N — l)bo is the length of the filament, yields the 
"activity number" A = Luq/k. The rates of active and 
elastic relaxation are = (Tq/'/jL'^ and = n/rjU^^^ 
respectively. Since A = Fc/Fk, the activity number also 
measures the ratio of time scales associated with active 
and elastic relaxation. As ^ — > the active time scale 
diverges and conformational changes occur only due to 
elastic forces. As ^ — >• oo conformational changes due to 
activity are much more rapid than those due to elasticity. 



Method: We use the lattice Boltzmann method [12] 
to obtain the incident flow in Eq. (3), corresponding to 
terms with m 7^ n, and then integrate the equations using 
the forward Euler method. The method of obtaining the 
incident flow is described in detail in [13] and in [14], 
and is related to methods used in [15]. The subtleties in 
using the lattice Boltzmann method to solve Eq. (3) are 
described in detail in [14]. We use lattice units in which 
both spatial and temporal discretization scales are unity. 
We choose 60 = 2 and k such that there is less than 1% 
variation in contour length. We choose R in the range 0.0 
to 0.5, ctq in the range to 0.05, and N in the range 16 to 
96. The initial filament conformation is a superposition 
of small amplitude transverse sinusoidal deformations of 
wavelengths a few integer multiples of L. The integration 
is carried out for several million time steps. Our results. 
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unless otherwise stated, are for periodic planar lattices 
of size 128. 

Results: We summarize our results in Figs. (1) and 
(2) and the movies in [14]. Our simulations find a 
linear instability of an initially straight filament. On 
dimensional grounds, this instability occurs only when 
L > Ia ^ k/cfq, where numerical prefactors can be ob- 
tained from a linear stability analysis of Eq. (3). In 
contrast, the elastic Euler instability of a filament un- 
der force F occurs when L > Ie ^ KjF . A linear 
instability of passive filaments in an active medium with- 
out nonlocal hydrodynamics was reported in [16], while 
bow-shaped conformation for filaments driven by exter- 
nal forces were reported in [17]. 

We show the nature of the linear instability, as .4 — > c« 
and fc — > 0, in Fig. (1) and its accompanying movie [14]. 
The flow produced by a straight filament is symmetric 
about the filament center and the filament axis, as shown 
in Fig. 1(a). It can only produce an extensile motion in 
the filament which is transient, since it is balanced by the 
stretching elasticity. A spontaneous transverse perturba- 
tion breaks the flow symmetry about the initial axis, en- 
hancing the perturbation, as shown in Fig. 1(b) . The flow 
now develops an uncompensated directional component 
in which the filament can translate. Since the hydrody- 
namic drag on the filament is greater at its ends [18], 
a balance between elastic deformation, active propulsion 
and drag ensues and the filament propels steadily in a 
deformed bow-shaped conformation [14]. 

In Fig. (1), the initial perturbation is symmetric about 
the filament midpoint. We call this an even mode. Initial 
perturbations which are antisymmetric about the mid- 
point are also linearly unstable. However, these odd 
modes produce flows of a completely different nature 
than the even modes. Instead of an uncompensated lin- 
ear component, the flow develops a vorticity centered on 
the filament midpoint which results in pure rotational 
motion [14]. A generic initial perturbation is a super- 
position of both even and odd modes and, thus, both 
rotates and translates. At low there is little cou- 
pling between the conformational modes, due to which 
the filament has steady motion. However, with increas- 
ing A greater elastic deformations are needed to balance 
the activity, due to which the filament develops nonlin- 
ear fiuctuating states with large- amplitude deformations, 
as shown in Fig. (2) and the accompanying movies [14]. 
Conformational modes are now coupled, and modes ab- 
sent from the initial perturbation can spontaneously ap- 
pear. With a predominance of even modes, the motion 
is primarily translational as seen in Fig. 2(a), but when 
a spontaneously generated odd mode appears the motion 
is both translational and rotational as seen in Fig. 2(b). 
In cubic flows, we see similar linear instabilities and non- 
linear fluctuating states [14]. It is surprising that such 
complex "animate" behavior can be generated by Eq. (3). 
Remarkably, its qualitative aspects can be understood 



from basic symmetry arguments relating conformations 
to the flows they produce. The role of symmetry in the 
motility of active droplets has been studied recently in 
[19]. 

In the free-draining approximation, it is possible to re- 
late the center of mass velocity R to the mean curvature 
of the fllament [14], 



where k is the local curvature and n is the local unit nor- 
mal vector. In Figs. 2(c) and 2(d) wc plot components of 
this equation corresponding to simulations in Figs. 2(a) 
and 2(b) respectively. The agreement is good in both 
cases. Thus local hydrodynamics provides a good ap- 
proximation for the translational velocity but non-local 
hydrodynamics is required to fully explain the confor- 
mational fluctuations. The interplay between nonlocal 
hydrodynamics and semiflexibility is necessary for the ro- 
tational motion of the fllament, as has been noted earlier 
in a different context [17]. 

For a microtubule of size L ^ 30fim, k BOpNum^ 
with about 200 motors per micron exerting approxi- 
mately 6pN of force, we obtain A ~ 60. These param- 
eters provide a translational speed of ~ lmms~^ for a 
semicircular shape in water. The activity can be manip- 
ulated in motor-polymer bundle systems or in polymers 
of autonomously motile nanorods over a large range of 
A. These systems would be the best candidates for a 
veriflcation of our results. 

Discussion and conclusion: Our model has several im- 
portant variations. We argued that active processes can- 
not add linear or angular momentum to the fluid and, 
so, must be represented by Stokesian singularities with 
those properties. This ruled out the Stokeslet and the 
rotlet but allowed for higher singularities, of which the 
stresslet, being the most dominant, was retained. The 
stresslet, with Sb Coo Q-xis, is not forbidden by symme- 
try as a representation of a polar active element. If it 
is forbidden for non-symmetry reasons, we must use the 
potential dipole d [10], the leading singularity with po- 
lar symmetry, whose velocity field is Va{r) = Gai3{'r)dp, 
Gap = (— (^Q^i +3faf^)/87r?7r^, in Eq. (3). The axis of the 
stresslet or the potential dipole can be oriented normally 
or obliquely to the local tangent of the filament and the 
stresslet can also be contractile, ctq < 0. The precise na- 
ture of the nonlinear steady states obtained from these 
various combination will be reported in future work. A 
generic equation of motion encompassing these specific 
cases is provided in [14]. 

Semiflexibility is crucially important in obtaining the 
results above. A rigid rod (k = oo, ^ = 0) is immune 
to the active instability. Since the uniaxial axes of the 
stresslets and the rod are aligned, it cannot, by symme- 
try, acquire any translational or rotational motion. It 
is only by the breaking of this symmetry, possible when 
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^ 7^ 0, that the filament is able to acquire motion. 

We have presented a model for an autonomously motile 
semiflexible filament which takes into account nonlocal 
hydrodynamic interactions. Our model opens up the pos- 
sibility of studying the nonequilibrium dynamics of active 
filaments, for example the cilia-like beating of motor- 
polymer bundles [7] as well as the collective properties 
of networks of active filaments, such as the cytoskeleton. 
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Supplemental material: 

Symmetries and flows : Spontaneous symmetry break- 
ing of the active filament generates net flows in which 
the filament can translate and/or rotate autonomously. 
In the absence of symmetry breaking, as seen in Fig. 
S3(a), the net flow on the filament cancels and no mo- 
tion is produced. Symmetry broken conformations can 
be distinguished as even or odd if they are symmetric 
or anti-symmetric, respectively, about an axis passing 
through the filament midpoint. In Fig. S3(b), an even 
conformation produces a net translatory fiow, with no 
rotational components. In Fig. S3(c), an odd conforma- 
tion produces a net rotatory fiow, with no translational 
components. The most general symmetry-broken confor- 
mation of the filament is a linear combination of even 
and odd conformations, allowing it to translate as well as 
rotate. 

Lattice Boltzmann method of solution : The lattice 
Boltzmann method (LBM) is an efficient way of solv- 
ing the Navier-Stokes equations. Here, we use it to solve 
for the fiow produced by the distribution of forces and 
stresses produced by the active filament in a periodic ge- 
ometry. The advantage of this is that it avoids the com- 
plicated Ewald summation of the hydrodynamic Green's 
functions. In particular, the LBM is easily extended 
to more complicated geometries, where closed form ex- 
pressions of the hydrodynamic Green's functions are not 
available. Specifically, we solve 

dtf + c, • V/, + [F • VJ], = -U, if, - /]>) (5) 

where fi (x, t) is the one-particle distribution function in 
phase space of coordinates x and velocities c^, Cij is the 
collision matrix linearized about the local equilibrium f^ 
and the repeated index j is summed over. The force 
density F is a sum of all the contributions f„ obtained 
from the potential energy of filament deformation and 
the stresslets on each monomer, modelled here as a pair 
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FIG. 3. FIG. SI. Effect of the spontaneous breaking of sym- 
metry on fluid flows, which for the motile (b) even and (c) 
odd symmetry broken conformations can be clearly distin- 
guished from that of (a), where symmetry is preserved and 
no autonomous motion is possible. Bead positions are shown 
as yellow filled circles, while the fluid velocity is shown as 
streamlines with the background color indicating its magni- 
tude. 
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of point forces separated by a distance d. This method- 
ology is explained in detail in [13]. The flow field at the 
location of the n-th bead, v(r„). is computed by inter- 
polation from the LBM grid points after "self" contri- 
butions from the monopole and dipole singularities have 
been subtracted. The bead is then moved forward by the 
equation of motion 



v(r„) 



(6) 



In this regard, our method is related to those of Diinweg 
and co-workers [15] and Fyta et al [15], where the LBM 
is used to solve for the fluid flow around a polymer. The 
main distinction is that those authors retain inertia for 
the beads while we work in the overdamped limit. 

To ensure that the LBM produces the Stokesian limit 
and that the velocity obtained above is identical to that 
produced by a direct summation of singularities, as in 
Eq. (3) of the main text, several conditions need to be 
ensured. First, the Reynolds number at the scale of the 
filament must be small, so that the nonlinear advection 
term in the Navier-Stokes equations is negligible. In our 
simulations, this is typically kept to about 10~^. Sec- 
ond, momentum diffusion across the filament must be 
rapid compared to the time scale of filament motion, so 
that there are no hydrodynamic retardation effects. In 
other words, the momentum diffusion time L'^/r] must be 
shorter than times of elastic and active relaxation. In our 
simulations, these ratios vary from 10~^ to 10~^. We al- 
ways ensure that the LBM operates in the hydrodynamic 
limit of small Knudsen number by choosing LBM relax- 
ation times (related to the shear and bulk viscosities) to 
be sufficiently below unity. 

Free-draining approximation of filament velocity : 
Consider the space curve, r — r{s). Taylor expanding 
about a point s gives, 

r(s') = r(s) + AsdsT + ^As'^d^r + 0{As^) (7) 

Using the Frenet-Scrret relations between the tangent 
t = dsT and the normal n we can recast the above ex- 
pansion as 

r{s') ^r{s) = Ast + ^As^>cn + 0{As^) (8) 

where the curvature x and torsion r are evaluated at s. 
If the curve is discretised by TV points, each separated 
by a distance bo, and inextensibility is enforced so that s 
remains a material parameter, we can write for the pairs 
of points r„_i, r„ and r„, r„+i 



1 2 

r„ - r„_i = 6ot„ - -b„xnn 
1 2 



(9) 
(10) 



The flows induced at r„ by stresslets at r„_i and r„+i, 
oriented along the local tangents, are ((To/47r?76Q)(r„ — 
r„_i) and (cro/47r?76Q)(r„ — r„_|_i) respectively. Utilising 
the above relations, the total flow at r„ due to neigh- 
bouring stresslets is 



r„, = - 



Ainjbo 



(11) 



In the nearest-neighbour approximation, the active ve- 
locity is proportional only to the curvature and is nor- 
mal to the curve. In the free-draining approximation, the 
equation of motion for an active fllament then is 



6Trrjbo imjbn 



(12) 



This can be thought of as the "Rouse" limit of Eq. (3) 
of the main text. The center of mass velocity can be ob- 
tained by summing all the bead velocities, while recalling 
that the sum of all forces is zero, to get 



R 



era 



inrjbo 



(13) 



Generic model for autonomously motile elastic fila- 
ments : We present a generic model for autonomously 
motile elastic fllaments which encompasses all the varia- 
tions of Eq. (3), discussed in the main text. We include 
the potential dipole as a possible singularity that is po- 
lar. This is subdominant to the stresslet, but is the most 
important singularity if the stresslet vanishes due to non- 
symmetry reasons. In addition, we include an externally 
imposed shear flow characterized by the shear rate ten- 
sor E. We allow for any orientation of the stresslet axis 
p„ and the potential dipole axis d„ relative to the local 
tangent of the filament t„. Thus, we have two new pa- 
rameters in the model, 9^ = p„ ■ t„ and 0^ = d„ • t„, the 
preferred angle that the stresslet and the potential dipole 
make with respect to the tangent. These angles can be 
made to vary along the filament, or may fluctuate in re- 
sponse to thermal noise. Finally, we include an external 
force g„ which may be due to externally imposed flelds 
like gravity or electricity. Such fields are required when 
studying the driven motion, for example sedimentation 
under gravity, of active filaments. The generic equation 
of motion then, is 
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N 



[0(r„ - r,„) ■ f,„ + D(r„ - r,„) • (t,„ + G(r„ - r,„) • d„] 



E 



(14) 



= - 



dU 
dr„ 



CTn = 0-o(PnPn - I/d), 



(15) 



0„/3(r) 



8m]r 



{Sap + fafp), Da/3yir) 



■{-faSpj + Sfafpfy), Gap = 



■{~Sap +ifafp) (16) 



For an unbounded two-dimensional fluid, the tensors can 
be obtained from their corresponding three-dimensional 
expressions through the replacements 1/r — >■ logr, 
l/r"+i 1/r" for n = 1,2 and 8tt An [9]. For pe- 
riodic flows, the forms given by Hasimoto must be used 
[20]. If required, the hydrodynamic interactions can be 
evaluated to higher orders in a multipolc expansion [21] 
or can be formulated within the more rigorous framework 
of slender body theory [18]. 

The relaxation rates associated with elasticity, stresslet 
activity and potential dipole activity are = k/t]L'^''^^, 
Fct = ao/riL'^ and F^ = do/rjL'^'^^. The shear rate tensor 
introduces at least one additional independent relaxation 
rate F^ = 7. The ratio of uniaxial and polar activities is 
Lao/ do, indicating that uniaxial activity dominates for 
long filaments. This motivates why Eq. (3) in the main 
text retains only uniaxial activity. These equations form 
the basis by which we can explore the noncquilibrium 
dynamics of active filaments under external fields or ex- 
ternally imposed velocity gradients. 
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